Intro

In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.

Load packages & source functions

# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(janitor)
library(forcats)
library(gapminder)
library(viridis)
library(ggridges)
library(raster)
library(icesDatras)
library(ggalluvial)
library(ggrepel)
library(ncdf4)
library(chron)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(geosphere)
library(quantreg)
#> Warning in .recacheSubclasses(def@className, def, env): undefined subclass
#> "numericVector" of class "Mnumeric"; definition not updated
library(brms)
#> Warning: package 'Rcpp' was built under R version 4.0.5
library(sdmTMB)
library(sp)
library(raster)
library(ggcorrplot)
library(ggh4x)
library(visreg)
options(mc.cores = parallel::detectCores()) 

# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/map_plot.R")

theme_set(theme_plot())

# Continuous colors
options(ggplot2.continuous.colour = "viridis")

# Discrete colors
scale_colour_discrete <- function(...) {
  scale_colour_brewer(palette = "Dark2")
}

scale_fill_discrete <- function(...) {
  scale_fill_brewer(palette = "Dark2")
}

Read and summarise data

cpue_full <- readr::read_csv("data/clean/catch_by_length_q1_q4.csv") %>%
  janitor::clean_names() %>% 
  rename(X = x,
         Y = y) 
#> Rows: 1759609 Columns: 22
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (8): Country, haul.id, IDx, ices_rect, id_haul_stomach, species, haul.i...
#> dbl (14): density, year, lat, lon, quarter, Month, sub_div, length_cm, depth...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)

# Summaries density by species group an haul
cpue <- cpue_full %>%
  mutate(species2 = species,
         species2 = ifelse(species == "cod" & length_cm >= 25, "large_cod", species),
         species2 = ifelse(species == "cod" & length_cm < 25, "small_cod", species2)) %>% 
  group_by(haul_id, species2) %>% 
  summarise(density = sum(density)) %>% 
  ungroup()
#> mutate: new variable 'species2' (character) with 3 unique values and 0% NA
#> group_by: 2 grouping variables (haul_id, species2)
#> summarise: now 27,897 rows and 3 columns, one group variable remaining (haul_id)
#> ungroup: no grouping variables

# Trim data with quantiles
ggplot(cpue, aes(density)) + 
  geom_histogram() +
  facet_wrap(~species2, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


cpue <- cpue %>% 
  group_by(species2) %>% 
  mutate(dens_quants = quantile(density, probs = 0.999)) %>% 
  filter(density < dens_quants) %>% 
  ungroup() %>% 
  dplyr::select(-dens_quants)
#> group_by: one grouping variable (species2)
#> mutate (grouped): new variable 'dens_quants' (double) with 3 unique values and 0% NA
#> filter (grouped): removed 30 rows (<1%), 27,867 rows remaining
#> ungroup: no grouping variables

# Make wide
wcpue <- cpue %>% pivot_wider(names_from = species2, values_from = density)
#> pivot_wider: reorganized (species2, density) into (flounder, large_cod, small_cod) [was 27867x3, now 9373x4]

head(wcpue)
#> # A tibble: 6 × 4
#>   haul_id                  flounder large_cod small_cod
#>   <chr>                       <dbl>     <dbl>     <dbl>
#> 1 1993:1:GFR:SOL:H20:21:1     9.41      77.9    0.628  
#> 2 1993:1:GFR:SOL:H20:22:32    6.63       5.13   0      
#> 3 1993:1:GFR:SOL:H20:23:31    0.953      2.61   0.00459
#> 4 1993:1:GFR:SOL:H20:24:30    1.89       9.71   0      
#> 5 1993:1:GFR:SOL:H20:25:2    16.7      400.     4.78   
#> 6 1993:1:GFR:SOL:H20:26:3    13.5      179.     2.95

# Left join in trawl information
nrow(wcpue)
#> [1] 9373
hauls <- cpue_full %>% distinct(haul_id, .keep_all = TRUE) %>% dplyr::select(-density)
#> distinct: removed 1,750,236 rows (99%), 9,373 rows remaining
nrow(hauls)
#> [1] 9373

cpue2 <- left_join(wcpue, hauls)
#> Joining, by = "haul_id"
#> left_join: added 20 columns (year, lat, lon, quarter, country, …)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 9,373
#> > =======
#> > rows total 9,373

colnames(cpue2)
#>  [1] "haul_id"         "flounder"        "large_cod"       "small_cod"      
#>  [5] "year"            "lat"             "lon"             "quarter"        
#>  [9] "country"         "month"           "i_dx"            "ices_rect"      
#> [13] "sub_div"         "length_cm"       "id_haul_stomach" "species"        
#> [17] "haul_id_size"    "substrate"       "depth"           "temp"           
#> [21] "oxy"             "sal"             "X"               "Y"

cpue2 <- cpue2 %>% drop_na(oxy, temp, depth, sal, substrate, flounder, small_cod, large_cod)
#> drop_na: removed 405 rows (4%), 8,968 rows remaining

# Scale variables
cpue2 <- cpue2 %>% 
  mutate(depth_ct = depth - mean(depth),
         depth_sq = depth_ct * depth_ct,
         depth_sq_sc = (depth_sq - mean(depth_sq)) / sd(depth_sq),
         depth_sc = (depth - mean(depth)) / sd(depth),
         temp_ct = temp - mean(temp),
         temp_sq = temp_ct * temp_ct,
         temp_sq_sc = (temp_sq - mean(temp_sq)) / sd(temp_sq),
         temp_sc = (temp - mean(temp)) / sd(temp),
         oxy_sc = (oxy - mean(oxy)) / sd(oxy),
         sal_sc = (sal - mean(sal)) / sd(sal),
         flounder_sc = (flounder - mean(flounder) / sd(flounder)),
         large_cod_sc = (large_cod - mean(large_cod) / sd(large_cod)),
         small_cod_sc = (small_cod - mean(small_cod) / sd(small_cod)),
         fyear = as.factor(year),
         fsubstrate = as.factor(substrate))
#> mutate: new variable 'depth_ct' (double) with 149 unique values and 0% NA
#>         new variable 'depth_sq' (double) with 149 unique values and 0% NA
#>         new variable 'depth_sq_sc' (double) with 149 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 149 unique values and 0% NA
#>         new variable 'temp_ct' (double) with 8,679 unique values and 0% NA
#>         new variable 'temp_sq' (double) with 8,679 unique values and 0% NA
#>         new variable 'temp_sq_sc' (double) with 8,679 unique values and 0% NA
#>         new variable 'temp_sc' (double) with 8,679 unique values and 0% NA
#>         new variable 'oxy_sc' (double) with 8,678 unique values and 0% NA
#>         new variable 'sal_sc' (double) with 8,690 unique values and 0% NA
#>         new variable 'flounder_sc' (double) with 7,933 unique values and 0% NA
#>         new variable 'large_cod_sc' (double) with 7,855 unique values and 0% NA
#>         new variable 'small_cod_sc' (double) with 6,617 unique values and 0% NA
#>         new variable 'fyear' (factor) with 28 unique values and 0% NA
#>         new variable 'fsubstrate' (factor) with 5 unique values and 0% NA
# colnames(cpue2)
# 
# ggplot(cpue2, aes(flounder, large_cod)) +
#   geom_point()
#   
# ggplot(cpue2, aes(flounder, small_cod)) +
#   geom_point()
# 
# write.csv(cpue2, "output/mich_catch_by_length_q1_q4.csv")

Read and scale the prediction grid

pred_grid_1 <- read_csv("data/clean/pred_grid_(1_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid_2 <- read_csv("data/clean/pred_grid_(2_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

pred_grid <- bind_rows(pred_grid_1, pred_grid_2)

# Scale with respect to data!
pred_grid <- pred_grid %>% 
  drop_na(substrate) %>% 
  mutate(X = X / 1000,
         Y = Y / 1000,
         depth_ct = depth - mean(cpue2$depth),
         depth_sq = depth_ct * depth_ct,
         depth_sq_sc = (depth_sq - mean(cpue2$depth_sq)) / sd(cpue2$depth_sq),
         depth_sc = (depth - mean(cpue2$depth)) / sd(cpue2$depth),
         temp_ct = temp - mean(cpue2$temp),
         temp_sq = temp_ct * temp_ct,
         temp_sq_sc = (temp_sq - mean(cpue2$temp_sq)) / sd(cpue2$temp_sq),
         temp_sc = (temp - mean(cpue2$temp)) / sd(cpue2$temp),
         oxy_sc = (oxy - mean(cpue2$oxy)) / sd(cpue2$oxy),
         sal_sc = (sal - mean(cpue2$sal)) / sd(cpue2$sal),
         fyear = as.factor(year),
         fsubstrate = as.factor(substrate))
#> drop_na: removed 280 rows (<1%), 592,760 rows remaining
#> mutate: changed 592,760 values (100%) of 'X' (0 new NA)
#>         changed 592,760 values (100%) of 'Y' (0 new NA)
#>         new variable 'depth_ct' (double) with 7,528 unique values and 0% NA
#>         new variable 'depth_sq' (double) with 7,528 unique values and 0% NA
#>         new variable 'depth_sq_sc' (double) with 7,528 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 7,528 unique values and 0% NA
#>         new variable 'temp_ct' (double) with 592,704 unique values and 0% NA
#>         new variable 'temp_sq' (double) with 592,704 unique values and 0% NA
#>         new variable 'temp_sq_sc' (double) with 592,704 unique values and 0% NA
#>         new variable 'temp_sc' (double) with 592,704 unique values and 0% NA
#>         new variable 'oxy_sc' (double) with 592,704 unique values and 0% NA
#>         new variable 'sal_sc' (double) with 592,760 unique values and 0% NA
#>         new variable 'fyear' (factor) with 28 unique values and 0% NA
#>         new variable 'fsubstrate' (factor) with 5 unique values and 0% NA

# Split by quarter
pred_grid_q1 <- pred_grid %>% filter(quarter == 1)
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
pred_grid_q4 <- pred_grid %>% filter(quarter == 4)
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
# Load models
# qwraps2::lazyload_cache_dir(path = "R/main_analysis/spatial_overlap_cache/html")

Plot data in space

# hist(cpue2$depth)
# hist(log(cpue2$depth))
# 
# cpue_long <- cpue2 %>%
#   pivot_longer(c(flounder, small_cod, large_cod)) %>% 
#   rename(species2 = name, density = value) %>% 
#   pivot_longer(c(depth, temp, oxy, sal)) %>% 
#   rename(env_var = name, env_var_value = value)
# 
# ggplot(cpue_long, aes(density)) + 
#   geom_histogram() + 
#   facet_wrap(~species2, scales = "free", ncol = 1)
# 
# ggplot(cpue_long, aes(env_var_value, density)) + 
#   geom_point(alpha = 0.3) + 
#   stat_smooth(method = "lm", formula = y ~ x, color = "blue", se = FALSE) + 
#   stat_smooth(method = "lm", formula = y ~ x + I(x^2), color = "red", se = FALSE) + 
#   facet_wrap(env_var~species2, scales = "free", ncol = 3) + 
#   coord_cartesian(expand = 0)
# 
# cpue2 <- cpue2 %>% mutate(fle_presence = ifelse(flounder == 0, "N", "Y"),
#                           l_cod_presence = ifelse(large_cod == 0, "N", "Y"),
#                           s_cod_presence = ifelse(small_cod == 0, "N", "Y"))
# 
# # Biomass density
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = flounder)) + 
#   scale_color_viridis_c(trans = "sqrt") + 
#   facet_wrap(~year)
# 
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = large_cod)) + 
#   scale_color_viridis_c(trans = "sqrt") + 
#   facet_wrap(~year)
# 
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = small_cod)) + 
#   scale_color_viridis_c(trans = "sqrt") + 
#   facet_wrap(~year)
# 
# # Presence / absence
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = fle_presence)) + 
#   facet_wrap(~year)
# 
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = l_cod_presence)) + 
#   facet_wrap(~year)
# 
# plot_map_fc + 
#   geom_point(data = cpue2, aes(X*1000, Y*1000, color = s_cod_presence)) + 
#   facet_wrap(~year)

Create meshes

# Both Q's
spde <- make_mesh(cpue2,
                  xy_cols = c("X", "Y"),
                  n_knots = 150,
                  type = "kmeans",
                  seed = 42)

# Q1
spde_q1 <- make_mesh(filter(cpue2, quarter == 1),
                     xy_cols = c("X", "Y"),
                     n_knots = 150, 
                     type = "kmeans",
                     seed = 42)
#> filter: removed 3,638 rows (41%), 5,330 rows remaining

# Q4
spde_q4 <- make_mesh(filter(cpue2, quarter == 4),
                     xy_cols = c("X", "Y"),
                     n_knots = 150, 
                     type = "kmeans",
                     seed = 42)
#> filter: removed 5,330 rows (59%), 3,638 rows remaining

Fit models

q1

# Small cod model q1 spatial
mcod_s_q1_space <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                            temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                          data = filter(cpue2, quarter == 1),
                          mesh = spde_q1,
                          family = tweedie(link = "log"),
                          spatiotemporal = "IID",
                          spatial = "on",
                          time = "year",
                          reml = FALSE,
                          control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,638 rows (41%), 5,330 rows remaining

sanity(mcod_s_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s_q1_space, conf.int = TRUE)
#>                             term    estimate  std.error   conf.low  conf.high
#> 1              fsubstratebedrock -0.54972097 1.08136018 -2.6691480  1.5697060
#> 2            fsubstratehard clay  0.17092271 0.67999722 -1.1618474  1.5036928
#> 3  fsubstratehard-bottom complex  0.03007446 0.70107772 -1.3440126  1.4041615
#> 4                  fsubstratemud -0.06888319 0.68011532 -1.4018847  1.2641183
#> 5                 fsubstratesand -0.18812057 0.68131988 -1.5234830  1.1472418
#> 6                      fyear1994 -0.27083562 0.65699125 -1.5585148  1.0168436
#> 7                      fyear1995  0.32913371 0.64556718 -0.9361547  1.5944221
#> 8                      fyear1996  0.20377555 0.65459642 -1.0792099  1.4867610
#> 9                      fyear1997 -0.78333814 0.65633408 -2.0697293  0.5030530
#> 10                     fyear1998 -0.15049380 0.64312981 -1.4110051  1.1100175
#> 11                     fyear1999 -0.34980464 0.63759208 -1.5994622  0.8998529
#> 12                     fyear2000  0.68322725 0.66498092 -0.6201114  1.9865659
#> 13                     fyear2001  0.53330300 0.63057208 -0.7025956  1.7692016
#> 14                     fyear2002  1.87881965 0.63487193  0.6344935  3.1231458
#> 15                     fyear2003  0.25683392 0.65287546 -1.0227785  1.5364463
#> 16                     fyear2004 -0.32113630 0.62569530 -1.5474766  0.9052040
#> 17                     fyear2005  1.94428972 0.60379069  0.7608817  3.1276977
#> 18                     fyear2006  0.17607793 0.61892330 -1.0369895  1.3891453
#> 19                     fyear2007  0.88916597 0.60511521 -0.2968381  2.0751700
#> 20                     fyear2008  1.22420903 0.60332687  0.0417101  2.4067080
#> 21                     fyear2009  0.74141683 0.60380274 -0.4420148  1.9248485
#> 22                     fyear2010  1.42101550 0.61003687  0.2253652  2.6166658
#> 23                     fyear2011  0.23229999 0.62086565 -0.9845743  1.4491743
#> 24                     fyear2012  0.33796827 0.61645874 -0.8702687  1.5462052
#> 25                     fyear2013  1.49256114 0.60443447  0.3078914  2.6772309
#> 26                     fyear2014  1.00617456 0.61271831 -0.1947313  2.2070804
#> 27                     fyear2015 -0.24874693 0.61446510 -1.4530764  0.9555825
#> 28                     fyear2016 -0.15828486 0.61471050 -1.3630953  1.0465256
#> 29                     fyear2017 -0.30106379 0.61950192 -1.5152652  0.9131377
#> 30                     fyear2018  0.64103439 0.61436656 -0.5631019  1.8451707
#> 31                     fyear2019 -0.15021962 0.66175918 -1.4472438  1.1468045
#> 32                     fyear2020 -1.20443537 0.66185914 -2.5016554  0.0927847
#> 33                      depth_sc  0.97037576 0.08787148  0.7981508  1.1426007
#> 34                   depth_sq_sc -1.28394502 0.07550417 -1.4319305 -1.1359596
#> 35                       temp_sc  1.01359298 0.13782015  0.7434705  1.2837155
#> 36                    temp_sq_sc -0.45935493 0.07285986 -0.6021576 -0.3165522
#> 37                        oxy_sc  0.58282738 0.08955571  0.4073014  0.7583533
#> 38                        sal_sc -0.08035481 0.10958924 -0.2951458  0.1344362
# Large cod model q1 spatial
mcod_l_q1_space <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                            temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                          data = filter(cpue2, quarter == 1),
                          mesh = spde_q1,
                          family = tweedie(link = "log"),
                          spatiotemporal = "IID",
                          spatial = "on",
                          time = "year",
                          reml = FALSE,
                          control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,638 rows (41%), 5,330 rows remaining

sanity(mcod_l_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l_q1_space, conf.int = TRUE)
#>                             term    estimate  std.error    conf.low   conf.high
#> 1              fsubstratebedrock  3.33190554 0.76118613  1.84000815  4.82380294
#> 2            fsubstratehard clay  3.96973799 0.47458439  3.03956967  4.89990631
#> 3  fsubstratehard-bottom complex  3.72985871 0.48918284  2.77107796  4.68863947
#> 4                  fsubstratemud  3.96832494 0.47482006  3.03769473  4.89895516
#> 5                 fsubstratesand  3.98132186 0.47566185  3.04904177  4.91360196
#> 6                      fyear1994  0.40602594 0.52096627 -0.61504919  1.42710107
#> 7                      fyear1995  0.12762982 0.51960316 -0.89077366  1.14603331
#> 8                      fyear1996  0.97962265 0.52204960 -0.04357576  2.00282106
#> 9                      fyear1997 -0.56368603 0.52591441 -1.59445932  0.46708727
#> 10                     fyear1998 -0.64602778 0.52346347 -1.67199734  0.37994177
#> 11                     fyear1999 -0.80355238 0.51875960 -1.82030251  0.21319775
#> 12                     fyear2000 -0.46274022 0.54747744 -1.53577629  0.61029585
#> 13                     fyear2001  0.04162165 0.51879413 -0.97519617  1.05843947
#> 14                     fyear2002  0.28583472 0.52835100 -0.74971420  1.32138365
#> 15                     fyear2003  0.26971748 0.53419605 -0.77728754  1.31672249
#> 16                     fyear2004 -1.53828619 0.51955950 -2.55660410 -0.51996827
#> 17                     fyear2005 -0.12409882 0.50031173 -1.10469180  0.85649416
#> 18                     fyear2006  0.35259729 0.50247871 -0.63224289  1.33743746
#> 19                     fyear2007 -0.28815905 0.50049926 -1.26911957  0.69280147
#> 20                     fyear2008  0.11722341 0.49765679 -0.85816597  1.09261280
#> 21                     fyear2009  0.38764681 0.49614749 -0.58478440  1.36007803
#> 22                     fyear2010  0.34810389 0.50198280 -0.63576431  1.33197210
#> 23                     fyear2011  0.09516022 0.50543903 -0.89548208  1.08580252
#> 24                     fyear2012 -0.45146797 0.50800862 -1.44714657  0.54421063
#> 25                     fyear2013 -0.13720362 0.49890510 -1.11503965  0.84063241
#> 26                     fyear2014 -0.23631283 0.50478973 -1.22568251  0.75305686
#> 27                     fyear2015 -0.16122928 0.50261088 -1.14632850  0.82386995
#> 28                     fyear2016 -0.22682356 0.50268367 -1.21206545  0.75841833
#> 29                     fyear2017 -0.54635722 0.50294671 -1.53211465  0.43940021
#> 30                     fyear2018 -0.38825368 0.50554395 -1.37910160  0.60259425
#> 31                     fyear2019 -1.02858635 0.54573504 -2.09820738  0.04103467
#> 32                     fyear2020 -1.30537254 0.53538228 -2.35470252 -0.25604256
#> 33                      depth_sc  0.75065191 0.06419549  0.62483106  0.87647276
#> 34                   depth_sq_sc -0.39245921 0.03479079 -0.46064789 -0.32427052
#> 35                       temp_sc  0.70706082 0.09830480  0.51438696  0.89973468
#> 36                    temp_sq_sc -0.14525828 0.05124781 -0.24570215 -0.04481441
#> 37                        oxy_sc  0.29743177 0.06620641  0.16766958  0.42719395
#> 38                        sal_sc -0.17440207 0.07928522 -0.32979825 -0.01900590
# Flounder model q1 spatial
mfle_q1_space <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                          temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                        data = filter(cpue2, quarter == 1),
                        mesh = spde_q1,
                        family = tweedie(link = "log"),
                        spatiotemporal = "IID",
                        spatial = "on",
                        time = "year",
                        reml = FALSE,
                        control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,638 rows (41%), 5,330 rows remaining

sanity(mfle_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle_q1_space, conf.int = TRUE)
#>                             term    estimate  std.error     conf.low
#> 1              fsubstratebedrock  4.09133618 0.52049350  3.071187672
#> 2            fsubstratehard clay  3.70825073 0.43066415  2.864164504
#> 3  fsubstratehard-bottom complex  3.60123785 0.43951943  2.739795601
#> 4                  fsubstratemud  3.78796875 0.43054351  2.944118973
#> 5                 fsubstratesand  3.68742673 0.43151917  2.841664697
#> 6                      fyear1994 -0.10648718 0.41370545 -0.917334955
#> 7                      fyear1995  0.03001146 0.41141847 -0.776353922
#> 8                      fyear1996  0.07201784 0.41606651 -0.743457541
#> 9                      fyear1997 -0.63738865 0.41839563 -1.457429017
#> 10                     fyear1998 -0.74615179 0.41391762 -1.557415421
#> 11                     fyear1999 -0.43988232 0.40810908 -1.239761423
#> 12                     fyear2000 -0.03650738 0.43029747 -0.879874913
#> 13                     fyear2001  0.07987207 0.41206573 -0.727761925
#> 14                     fyear2002  0.37427374 0.41743409 -0.443882046
#> 15                     fyear2003  0.28430880 0.41847879 -0.535894551
#> 16                     fyear2004 -0.95578850 0.40545811 -1.750471793
#> 17                     fyear2005  0.05309166 0.39438861 -0.719895802
#> 18                     fyear2006  0.60071762 0.39551691 -0.174481285
#> 19                     fyear2007  0.34827399 0.39233210 -0.420682810
#> 20                     fyear2008  0.27884936 0.39317426 -0.491758030
#> 21                     fyear2009  0.14276616 0.39269484 -0.626901588
#> 22                     fyear2010  0.66332645 0.39474908 -0.110367540
#> 23                     fyear2011  0.21630668 0.39770089 -0.563172740
#> 24                     fyear2012  0.10076256 0.39578708 -0.674965867
#> 25                     fyear2013  0.41927358 0.39221728 -0.349458156
#> 26                     fyear2014  0.19635692 0.39694484 -0.581640683
#> 27                     fyear2015 -0.27416323 0.39674972 -1.051778392
#> 28                     fyear2016 -0.12358694 0.39610329 -0.899935112
#> 29                     fyear2017 -0.08079854 0.39553248 -0.856027961
#> 30                     fyear2018  0.26841628 0.39925996 -0.514118857
#> 31                     fyear2019 -0.14305015 0.42921210 -0.984290404
#> 32                     fyear2020 -0.29538565 0.41996741 -1.118506654
#> 33                      depth_sc  0.09669746 0.05062436 -0.002524471
#> 34                   depth_sq_sc -0.09632404 0.02637195 -0.148012125
#> 35                       temp_sc  0.81118129 0.08580203  0.643012399
#> 36                    temp_sq_sc -0.02569753 0.04448537 -0.112887244
#> 37                        oxy_sc  0.51694645 0.05389546  0.411313283
#> 38                        sal_sc  0.25730289 0.07541040  0.109501221
#>      conf.high
#> 1   5.11148469
#> 2   4.55233697
#> 3   4.46268010
#> 4   4.63181852
#> 5   4.53318877
#> 6   0.70436060
#> 7   0.83637685
#> 8   0.88749322
#> 9   0.18265172
#> 10  0.06511185
#> 11  0.35999679
#> 12  0.80686016
#> 13  0.88750606
#> 14  1.19242952
#> 15  1.10451215
#> 16 -0.16110521
#> 17  0.82607913
#> 18  1.37591652
#> 19  1.11723078
#> 20  1.04945675
#> 21  0.91243391
#> 22  1.43702043
#> 23  0.99578610
#> 24  0.87649099
#> 25  1.18800532
#> 26  0.97435451
#> 27  0.50345194
#> 28  0.65276124
#> 29  0.69443088
#> 30  1.05095142
#> 31  0.69819011
#> 32  0.52773536
#> 33  0.19591939
#> 34 -0.04463596
#> 35  0.97935018
#> 36  0.06149219
#> 37  0.62257961
#> 38  0.40510457

q4

# Small cod model q4 spatial
mcod_s_q4_space <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                            temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                          data = filter(cpue2, quarter == 4),
                          mesh = spde_q4,
                          family = tweedie(link = "log"),
                          spatiotemporal = "IID",
                          spatial = "on",
                          time = "year",
                          reml = FALSE,
                          control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,330 rows (59%), 3,638 rows remaining

sanity(mcod_s_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s_q4_space, conf.int = TRUE)
#>                             term    estimate  std.error   conf.low   conf.high
#> 1              fsubstratebedrock  0.32249199 0.73538955 -1.1188450  1.76382903
#> 2            fsubstratehard clay -0.62378219 0.55544819 -1.7124406  0.46487626
#> 3  fsubstratehard-bottom complex -0.96634332 0.57126621 -2.0860045  0.15331788
#> 4                  fsubstratemud -0.89407109 0.55349970 -1.9789106  0.19076839
#> 5                 fsubstratesand -0.67373306 0.55549607 -1.7624853  0.41501922
#> 6                      fyear1994  0.06757994 0.68460289 -1.2742171  1.40937695
#> 7                      fyear1995  1.05387153 0.65331032 -0.2265932  2.33433623
#> 8                      fyear1996 -0.74286053 0.69564320 -2.1062961  0.62057508
#> 9                      fyear1997  0.54977163 0.63147742 -0.6879014  1.78744463
#> 10                     fyear1998  0.31989404 0.65061515 -0.9552882  1.59507631
#> 11                     fyear1999  0.51812992 0.63184376 -0.7202611  1.75652093
#> 12                     fyear2000  1.07924843 0.60720418 -0.1108499  2.26934676
#> 13                     fyear2001  2.04473353 0.59447261  0.8795886  3.20987843
#> 14                     fyear2002  1.44874905 0.59584509  0.2809141  2.61658397
#> 15                     fyear2003 -0.30863035 0.61053166 -1.5052504  0.88798971
#> 16                     fyear2004  2.37016410 0.58870662  1.2163203  3.52400788
#> 17                     fyear2005  1.47212495 0.57301842  0.3490295  2.59522040
#> 18                     fyear2006  1.83767110 0.56994858  0.7205924  2.95474978
#> 19                     fyear2007  2.00583510 0.56804833  0.8924808  3.11918938
#> 20                     fyear2008  1.54410451 0.56548803  0.4357683  2.65244068
#> 21                     fyear2009  2.26931005 0.56514325  1.1616496  3.37697046
#> 22                     fyear2010  1.57705454 0.56759490  0.4645890  2.68952010
#> 23                     fyear2011  0.75334444 0.57103716 -0.3658678  1.87255671
#> 24                     fyear2012  1.31095911 0.57701625  0.1800280  2.44189018
#> 25                     fyear2013  1.86685633 0.57015276  0.7493775  2.98433522
#> 26                     fyear2014  0.73565198 0.57489336 -0.3911183  1.86242226
#> 27                     fyear2015  0.33706211 0.57855532 -0.7968855  1.47100970
#> 28                     fyear2016 -0.01423134 0.57677617 -1.1446919  1.11622918
#> 29                     fyear2017  1.41721173 0.56912509  0.3017471  2.53267640
#> 30                     fyear2018  0.37104157 0.57473009 -0.7554087  1.49749184
#> 31                     fyear2019 -0.15607354 0.63784335 -1.4062235  1.09407646
#> 32                     fyear2020  0.12153068 0.60740905 -1.0689692  1.31203054
#> 33                      depth_sc  0.28276560 0.09106935  0.1042730  0.46125823
#> 34                   depth_sq_sc -1.36146413 0.08572625 -1.5294845 -1.19344378
#> 35                       temp_sc  0.82467933 0.13390142  0.5622374  1.08712129
#> 36                    temp_sq_sc -0.20336686 0.07336668 -0.3471629 -0.05957081
#> 37                        oxy_sc  0.99602212 0.10423533  0.7917246  1.20031962
#> 38                        sal_sc  0.46967220 0.11623198  0.2418617  0.69748268
# Large cod model q4 spatial
mcod_l_q4_space <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                            temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                          data = filter(cpue2, quarter == 4),
                          mesh = spde_q4,
                          family = tweedie(link = "log"),
                          spatiotemporal = "IID",
                          spatial = "on",
                          time = "year",
                          reml = FALSE,
                          control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,330 rows (59%), 3,638 rows remaining

sanity(mcod_l_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l_q4_space, conf.int = TRUE)
#>                             term     estimate  std.error    conf.low
#> 1              fsubstratebedrock  4.461261842 0.67727269  3.13383176
#> 2            fsubstratehard clay  4.085728133 0.50076507  3.10424663
#> 3  fsubstratehard-bottom complex  3.766673734 0.51197346  2.76322418
#> 4                  fsubstratemud  3.913855555 0.49969560  2.93447018
#> 5                 fsubstratesand  3.887435320 0.50143270  2.90464528
#> 6                      fyear1994  0.419584087 0.56203685 -0.68198790
#> 7                      fyear1995  1.069600130 0.55806448 -0.02418614
#> 8                      fyear1996  0.257708655 0.56677971 -0.85315916
#> 9                      fyear1997 -0.370530757 0.54713854 -1.44290259
#> 10                     fyear1998 -0.204480195 0.55538387 -1.29301257
#> 11                     fyear1999 -0.631959902 0.55349467 -1.71678952
#> 12                     fyear2000  0.073268003 0.53256216 -0.97053465
#> 13                     fyear2001  0.160909371 0.52872110 -0.87536495
#> 14                     fyear2002  0.776504998 0.52166877 -0.24594700
#> 15                     fyear2003 -1.483437852 0.53590092 -2.53378435
#> 16                     fyear2004  0.096387087 0.52463275 -0.93187420
#> 17                     fyear2005  0.600156755 0.50644354 -0.39245434
#> 18                     fyear2006  0.332251050 0.50594943 -0.65939161
#> 19                     fyear2007  0.512317962 0.50465526 -0.47678817
#> 20                     fyear2008  0.655240592 0.50146897 -0.32762053
#> 21                     fyear2009  0.513073198 0.50309145 -0.47296793
#> 22                     fyear2010  0.676087205 0.50348831 -0.31073176
#> 23                     fyear2011 -0.007275424 0.50555201 -0.99813915
#> 24                     fyear2012 -0.396814061 0.51432782 -1.40487806
#> 25                     fyear2013  0.002090184 0.50732834 -0.99225509
#> 26                     fyear2014 -0.202627548 0.50870179 -1.19966474
#> 27                     fyear2015 -0.087591060 0.50823553 -1.08371439
#> 28                     fyear2016 -0.818491405 0.51022386 -1.81851179
#> 29                     fyear2017 -0.424939151 0.50767462 -1.41996313
#> 30                     fyear2018 -1.153193343 0.51144123 -2.15559973
#> 31                     fyear2019 -1.254966349 0.55782243 -2.34827821
#> 32                     fyear2020 -1.536938007 0.54214527 -2.59952322
#> 33                      depth_sc  0.408087994 0.07078593  0.26935012
#> 34                   depth_sq_sc -0.542902793 0.05660972 -0.65385581
#> 35                       temp_sc  0.239683585 0.10150667  0.04073417
#> 36                    temp_sq_sc -0.063530409 0.05595097 -0.17319229
#> 37                        oxy_sc  0.722704596 0.07857480  0.56870083
#> 38                        sal_sc  0.288968482 0.09016190  0.11225441
#>      conf.high
#> 1   5.78869192
#> 2   5.06720964
#> 3   4.77012329
#> 4   4.89324093
#> 5   4.87022536
#> 6   1.52115608
#> 7   2.16338640
#> 8   1.36857647
#> 9   0.70184108
#> 10  0.88405218
#> 11  0.45286972
#> 12  1.11707066
#> 13  1.19718369
#> 14  1.79895699
#> 15 -0.43309135
#> 16  1.12464838
#> 17  1.59276785
#> 18  1.32389371
#> 19  1.50142409
#> 20  1.63810171
#> 21  1.49911433
#> 22  1.66290617
#> 23  0.98358831
#> 24  0.61124994
#> 25  0.99643546
#> 26  0.79440965
#> 27  0.90853227
#> 28  0.18152898
#> 29  0.57008483
#> 30 -0.15078695
#> 31 -0.16165448
#> 32 -0.47435279
#> 33  0.54682587
#> 34 -0.43194978
#> 35  0.43863300
#> 36  0.04613147
#> 37  0.87670836
#> 38  0.46568256
# Flounder model q4 spatial
mfle_q4_space <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                          temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                      data = filter(cpue2, quarter == 4),
                      mesh = spde_q4,
                      family = tweedie(link = "log"),
                      spatiotemporal = "IID",
                      spatial = "on",
                      time = "year",
                      reml = FALSE,
                      control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,330 rows (59%), 3,638 rows remaining

sanity(mfle_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle_q4_space, conf.int = TRUE)
#>                             term    estimate  std.error    conf.low   conf.high
#> 1              fsubstratebedrock  2.71204220 0.60223670  1.53167996  3.89240444
#> 2            fsubstratehard clay  2.09065120 0.56330160  0.98660035  3.19470205
#> 3  fsubstratehard-bottom complex  1.98083675 0.56839418  0.86680462  3.09486888
#> 4                  fsubstratemud  2.01546063 0.56224288  0.91348483  3.11743643
#> 5                 fsubstratesand  2.10477305 0.56324257  1.00083789  3.20870820
#> 6                      fyear1994 -0.42238110 0.59608389 -1.59068405  0.74592184
#> 7                      fyear1995  0.82691107 0.58064409 -0.31113043  1.96495257
#> 8                      fyear1996  0.34292251 0.59481583 -0.82289509  1.50874011
#> 9                      fyear1997  0.24815800 0.56532853 -0.85986556  1.35618156
#> 10                     fyear1998 -0.26640178 0.57547313 -1.39430839  0.86150484
#> 11                     fyear1999 -0.96023820 0.57740366 -2.09192858  0.17145219
#> 12                     fyear2000  0.06955067 0.55504887 -1.01832512  1.15742646
#> 13                     fyear2001  0.14979012 0.55147291 -0.93107693  1.23065717
#> 14                     fyear2002  0.87163543 0.54406255 -0.19470757  1.93797844
#> 15                     fyear2003 -0.99038344 0.54565467 -2.05984694  0.07908007
#> 16                     fyear2004  0.51169827 0.54311299 -0.55278363  1.57618016
#> 17                     fyear2005  0.52266141 0.52661866 -0.50949221  1.55481502
#> 18                     fyear2006  0.58363248 0.52306129 -0.44154881  1.60881377
#> 19                     fyear2007  0.68479501 0.52247358 -0.33923438  1.70882440
#> 20                     fyear2008  0.88301325 0.51990783 -0.13598738  1.90201388
#> 21                     fyear2009  0.90246533 0.51882379 -0.11441062  1.91934127
#> 22                     fyear2010  1.21271007 0.51853304  0.19640400  2.22901615
#> 23                     fyear2011  0.87184325 0.51905899 -0.14549368  1.88918019
#> 24                     fyear2012  0.48386735 0.52304030 -0.54127280  1.50900751
#> 25                     fyear2013  0.71014527 0.52389239 -0.31666494  1.73695548
#> 26                     fyear2014  0.42636457 0.52325026 -0.59918710  1.45191624
#> 27                     fyear2015  0.49211851 0.52179361 -0.53057818  1.51481520
#> 28                     fyear2016  0.36022963 0.52073115 -0.66038467  1.38084393
#> 29                     fyear2017  0.56235857 0.52240042 -0.46152744  1.58624458
#> 30                     fyear2018  0.37786922 0.52380208 -0.64876400  1.40450244
#> 31                     fyear2019  0.58021712 0.56227654 -0.52182464  1.68225888
#> 32                     fyear2020  0.03493564 0.54492131 -1.03309051  1.10296178
#> 33                      depth_sc -0.06312977 0.06548158 -0.19147132  0.06521178
#> 34                   depth_sq_sc -0.83253692 0.06332174 -0.95664524 -0.70842859
#> 35                       temp_sc  0.57487192 0.10083785  0.37723337  0.77251048
#> 36                    temp_sq_sc -0.03971065 0.05665936 -0.15076096  0.07133966
#> 37                        oxy_sc  1.20536646 0.08323146  1.04223581  1.36849712
#> 38                        sal_sc  0.16223809 0.09127951 -0.01666646  0.34114264

Models with the other species as coeffcients

We can’t put them in the main model because we want to use that for spatial overlap. There isn’t as much need to fit the models below by quarter (though we could), because while the distribution changes, we don’t know a priori that the effect should.

# Small cod model
mcod_s <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                         temp_sc + temp_sq_sc + oxy_sc + sal_sc + flounder_sc,
                       data = cpue2,
                       mesh = spde,
                       family = tweedie(link = "log"),
                       spatiotemporal = "IID",
                       spatial = "on",
                       time = "year",
                       reml = FALSE,
                       control = sdmTMBcontrol(newton_loops = 1))
#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation

#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation

sanity(mcod_s)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s, conf.int = TRUE)
#>                             term     estimate    std.error    conf.low
#> 1              fsubstratebedrock  0.628988603 0.6284130899 -0.60267842
#> 2            fsubstratehard clay  0.186705195 0.4303576376 -0.65678028
#> 3  fsubstratehard-bottom complex  0.139329862 0.4463497732 -0.73549962
#> 4                  fsubstratemud  0.021906096 0.4298307553 -0.82054670
#> 5                 fsubstratesand  0.045164068 0.4307133542 -0.79901859
#> 6                      fyear1994 -0.647472621 0.4659866755 -1.56078972
#> 7                      fyear1995  0.135138774 0.4538600093 -0.75441050
#> 8                      fyear1996 -0.585039168 0.4663404873 -1.49904973
#> 9                      fyear1997 -0.781861845 0.4576482420 -1.67883592
#> 10                     fyear1998 -0.456436267 0.4538632659 -1.34599192
#> 11                     fyear1999 -0.475663483 0.4499229927 -1.35749634
#> 12                     fyear2000  0.363354134 0.4566324637 -0.53162905
#> 13                     fyear2001  0.581528036 0.4369675478 -0.27491262
#> 14                     fyear2002  1.002107172 0.4341469713  0.15119474
#> 15                     fyear2003 -0.503903733 0.4517292266 -1.38927675
#> 16                     fyear2004  0.432451462 0.4312379975 -0.41275948
#> 17                     fyear2005  1.149304640 0.4156855586  0.33457592
#> 18                     fyear2006  0.496325702 0.4223847818 -0.33153326
#> 19                     fyear2007  0.759332204 0.4169238409 -0.05782351
#> 20                     fyear2008  0.771594631 0.4156141740 -0.04299418
#> 21                     fyear2009  0.990532701 0.4140889987  0.17893318
#> 22                     fyear2010  0.831369822 0.4184572443  0.01120869
#> 23                     fyear2011 -0.040597968 0.4221376512 -0.86797256
#> 24                     fyear2012  0.251791028 0.4241668674 -0.57956076
#> 25                     fyear2013  1.225870027 0.4168090620  0.40893928
#> 26                     fyear2014  0.546055609 0.4232892680 -0.28357611
#> 27                     fyear2015 -0.426740925 0.4262169679 -1.26211083
#> 28                     fyear2016 -0.370330338 0.4232011663 -1.19978938
#> 29                     fyear2017  0.309418397 0.4208926310 -0.51551600
#> 30                     fyear2018  0.135837741 0.4199787271 -0.68730544
#> 31                     fyear2019 -0.404787529 0.4623097601 -1.31089801
#> 32                     fyear2020 -1.015116571 0.4529117231 -1.90280724
#> 33                      depth_sc  0.563637368 0.0671644517  0.43199746
#> 34                   depth_sq_sc -1.086366763 0.0591121348 -1.20222442
#> 35                       temp_sc  0.620151439 0.0304481960  0.56047407
#> 36                    temp_sq_sc -0.380384042 0.0230597138 -0.42558025
#> 37                        oxy_sc  0.506692384 0.0572793122  0.39442700
#> 38                        sal_sc  0.070215044 0.0591226458 -0.04566321
#> 39                   flounder_sc  0.002697162 0.0001541721  0.00239499
#>       conf.high
#> 1   1.860655627
#> 2   1.030190665
#> 3   1.014159342
#> 4   0.864358896
#> 5   0.889346730
#> 6   0.265844480
#> 7   1.024688046
#> 8   0.328971392
#> 9   0.115112227
#> 10  0.433119388
#> 11  0.406169378
#> 12  1.258337317
#> 13  1.437968692
#> 14  1.853019599
#> 15  0.381469281
#> 16  1.277662405
#> 17  1.964033363
#> 18  1.324184662
#> 19  1.576487917
#> 20  1.586183444
#> 21  1.802132225
#> 22  1.651530950
#> 23  0.786776625
#> 24  1.083142811
#> 25  2.042800777
#> 26  1.375687329
#> 27  0.408628982
#> 28  0.459128706
#> 29  1.134352795
#> 30  0.958980920
#> 31  0.501322951
#> 32 -0.127425906
#> 33  0.695277274
#> 34 -0.970509107
#> 35  0.679828807
#> 36 -0.335187834
#> 37  0.618957773
#> 38  0.186093301
#> 39  0.002999333
# Large cod model q4 spatial
mcod_l <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                   temp_sc + temp_sq_sc + oxy_sc + sal_sc + flounder_sc,
                 data = cpue2,
                 mesh = spde,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 spatial = "on",
                 time = "year",
                 reml = FALSE,
                 control = sdmTMBcontrol(newton_loops = 1))
#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation

sanity(mcod_l)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l, conf.int = TRUE)
#>                             term    estimate    std.error     conf.low
#> 1              fsubstratebedrock  3.42710819 0.5776955169  2.294845779
#> 2            fsubstratehard clay  3.65473107 0.4329649614  2.806135343
#> 3  fsubstratehard-bottom complex  3.53484686 0.4407693938  2.670954722
#> 4                  fsubstratemud  3.63725350 0.4326850367  2.789206409
#> 5                 fsubstratesand  3.61916586 0.4331881805  2.770132626
#> 6                      fyear1994  0.19742496 0.4461982359 -0.677107517
#> 7                      fyear1995  0.17722202 0.4452049165 -0.695363579
#> 8                      fyear1996  0.59342116 0.4467822919 -0.282256041
#> 9                      fyear1997 -0.62579353 0.4465154309 -1.500947688
#> 10                     fyear1998 -0.58927473 0.4457547961 -1.462938074
#> 11                     fyear1999 -0.78880418 0.4439619893 -1.658953688
#> 12                     fyear2000 -0.16735879 0.4541355806 -1.057448169
#> 13                     fyear2001 -0.10073693 0.4411441110 -0.965363495
#> 14                     fyear2002  0.29062836 0.4386040328 -0.569019749
#> 15                     fyear2003 -0.59834008 0.4507702144 -1.481833463
#> 16                     fyear2004 -0.68491754 0.4355127979 -1.538506941
#> 17                     fyear2005  0.16123839 0.4238126275 -0.669419100
#> 18                     fyear2006  0.30196408 0.4258227266 -0.532633127
#> 19                     fyear2007 -0.05063772 0.4254537613 -0.884511765
#> 20                     fyear2008  0.23752649 0.4236463810 -0.592805160
#> 21                     fyear2009  0.47468725 0.4219802066 -0.352378757
#> 22                     fyear2010  0.47630118 0.4249369173 -0.356559877
#> 23                     fyear2011 -0.10742847 0.4265686555 -0.943487677
#> 24                     fyear2012 -0.45734127 0.4305572091 -1.301217898
#> 25                     fyear2013 -0.06534330 0.4252864066 -0.898889341
#> 26                     fyear2014 -0.10135916 0.4285620124 -0.941325273
#> 27                     fyear2015 -0.18794203 0.4278915386 -1.026594032
#> 28                     fyear2016 -0.55035878 0.4276565701 -1.388550257
#> 29                     fyear2017 -0.52852784 0.4268257001 -1.365090842
#> 30                     fyear2018 -0.72796029 0.4276077623 -1.566056107
#> 31                     fyear2019 -1.15141866 0.4604934626 -2.053969258
#> 32                     fyear2020 -1.47517542 0.4485627645 -2.354342285
#> 33                      depth_sc  0.48942278 0.0492980160  0.392800447
#> 34                   depth_sq_sc -0.31765342 0.0284809234 -0.373475003
#> 35                       temp_sc  0.32629619 0.0213972972  0.284358257
#> 36                    temp_sq_sc -0.25106483 0.0164175004 -0.283242543
#> 37                        oxy_sc  0.23252038 0.0414262955  0.151326338
#> 38                        sal_sc -0.14439376 0.0426861689 -0.228057118
#> 39                   flounder_sc  0.00223198 0.0001122465  0.002011981
#>       conf.high
#> 1   4.559370594
#> 2   4.503326805
#> 3   4.398738997
#> 4   4.485300587
#> 5   4.468199091
#> 6   1.071957428
#> 7   1.049807625
#> 8   1.469098361
#> 9   0.249360638
#> 10  0.284388618
#> 11  0.081345331
#> 12  0.722730595
#> 13  0.763889644
#> 14  1.150276466
#> 15  0.285153308
#> 16  0.168671857
#> 17  0.991895872
#> 18  1.136561289
#> 19  0.783236333
#> 20  1.067858138
#> 21  1.301753257
#> 22  1.309162231
#> 23  0.728630727
#> 24  0.386535348
#> 25  0.768202739
#> 26  0.738606946
#> 27  0.650709978
#> 28  0.287832693
#> 29  0.308035157
#> 30  0.110135520
#> 31 -0.248868055
#> 32 -0.596008558
#> 33  0.586045118
#> 34 -0.261831835
#> 35  0.368234121
#> 36 -0.218887124
#> 37  0.313714432
#> 38 -0.060730410
#> 39  0.002451979
# Flounder model q4 spatial
mfle <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                 temp_sc + temp_sq_sc + oxy_sc + sal_sc + large_cod_sc + small_cod_sc,
               data = cpue2,
               mesh = spde,
               family = tweedie(link = "log"),
               spatiotemporal = "IID",
               spatial = "on",
               time = "year",
               reml = FALSE,
               control = sdmTMBcontrol(newton_loops = 1))
#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation

#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation

sanity(mfle)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle, conf.int = TRUE)
#>                             term      estimate    std.error      conf.low
#> 1              fsubstratebedrock  3.7480605290 4.809351e-01  2.8054450987
#> 2            fsubstratehard clay  3.3821101833 4.265386e-01  2.5461098645
#> 3  fsubstratehard-bottom complex  3.4792137737 4.328010e-01  2.6309394082
#> 4                  fsubstratemud  3.5188278644 4.260961e-01  2.6836949097
#> 5                 fsubstratesand  3.4790847848 4.272366e-01  2.6417164564
#> 6                      fyear1994 -0.1746850881 3.764376e-01 -0.9124892955
#> 7                      fyear1995  0.0322147643 3.746166e-01 -0.7020202906
#> 8                      fyear1996  0.0789515893 3.779533e-01 -0.6618232722
#> 9                      fyear1997 -0.5308645488 3.770456e-01 -1.2698602707
#> 10                     fyear1998 -0.7059549558 3.762681e-01 -1.4434268281
#> 11                     fyear1999 -0.5208658197 3.709189e-01 -1.2478534311
#> 12                     fyear2000 -0.2732915952 3.851104e-01 -1.0280940453
#> 13                     fyear2001 -0.0062696627 3.727832e-01 -0.7369112902
#> 14                     fyear2002  0.2945432822 3.674595e-01 -0.4256641871
#> 15                     fyear2003 -0.3504426411 3.734771e-01 -1.0824442460
#> 16                     fyear2004 -0.3602718594 3.635209e-01 -1.0727596604
#> 17                     fyear2005 -0.0703193795 3.544560e-01 -0.7650404426
#> 18                     fyear2006  0.3092207848 3.551067e-01 -0.3867755284
#> 19                     fyear2007  0.1500276245 3.534010e-01 -0.5426256724
#> 20                     fyear2008  0.2362681297 3.535750e-01 -0.4567262252
#> 21                     fyear2009  0.1724025532 3.531042e-01 -0.5196688931
#> 22                     fyear2010  0.5688610670 3.539224e-01 -0.1248140049
#> 23                     fyear2011  0.1819444743 3.544701e-01 -0.5128041571
#> 24                     fyear2012  0.0013916538 3.562480e-01 -0.6968415743
#> 25                     fyear2013  0.3524349089 3.537225e-01 -0.3408483580
#> 26                     fyear2014 -0.1340147825 3.581498e-01 -0.8359754149
#> 27                     fyear2015 -0.2295553924 3.567940e-01 -0.9288587597
#> 28                     fyear2016 -0.2423025107 3.555494e-01 -0.9391664435
#> 29                     fyear2017 -0.1301990053 3.548341e-01 -0.8256611383
#> 30                     fyear2018  0.0039169536 3.565250e-01 -0.6948592487
#> 31                     fyear2019  0.0590546671 3.854075e-01 -0.6963301850
#> 32                     fyear2020 -0.4085526520 3.748436e-01 -1.1432325475
#> 33                      depth_sc  0.1779467626 4.925036e-02  0.0814178360
#> 34                   depth_sq_sc -0.0236437979 2.638278e-02 -0.0753531027
#> 35                       temp_sc  0.4247440518 2.082157e-02  0.3839345247
#> 36                    temp_sq_sc -0.3098803388 1.650299e-02 -0.3422256124
#> 37                        oxy_sc  0.3911327055 4.277491e-02  0.3072954156
#> 38                        sal_sc  0.2797308385 4.675097e-02  0.1881006181
#> 39                  large_cod_sc  0.0005973076 4.037953e-05  0.0005181651
#> 40                  small_cod_sc  0.0038476870 4.201599e-04  0.0030241887
#>       conf.high
#> 1   4.690675959
#> 2   4.218110502
#> 3   4.327488139
#> 4   4.353960819
#> 5   4.316453113
#> 6   0.563119119
#> 7   0.766449819
#> 8   0.819726451
#> 9   0.208131173
#> 10  0.031516917
#> 11  0.206121792
#> 12  0.481510855
#> 13  0.724371965
#> 14  1.014750751
#> 15  0.381558964
#> 16  0.352215942
#> 17  0.624401684
#> 18  1.005217098
#> 19  0.842680922
#> 20  0.929262485
#> 21  0.864474000
#> 22  1.262536139
#> 23  0.876693106
#> 24  0.699624882
#> 25  1.045718176
#> 26  0.567945850
#> 27  0.469747975
#> 28  0.454561422
#> 29  0.565263128
#> 30  0.702693156
#> 31  0.814439519
#> 32  0.326127243
#> 33  0.274475689
#> 34  0.028065507
#> 35  0.465553579
#> 36 -0.277535065
#> 37  0.474969995
#> 38  0.371361059
#> 39  0.000676450
#> 40  0.004671185

Plot conditional density effects

mfle_d <- visreg(mfle, xvar = "small_cod_sc", plot = FALSE)
mfle_d2 <- visreg(mfle, xvar = "large_cod_sc", plot = FALSE)

mscod_d <- visreg(mcod_s, xvar = "flounder_sc", plot = FALSE)
mlcod_d <- visreg(mcod_l, xvar = "flounder_sc", plot = FALSE)
pal <- brewer.pal(n = 3, name = "Set1")[2]

p1 <- ggplot(mfle_d$fit, aes(x = small_cod_sc, y = visregFit)) +
  geom_point(aes(y = visregRes), data = mfle_d$res, size = 1, alpha = 0.2, color = "grey30") +
  geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2) +
  geom_line(color = pal, size = 1.2) +
  labs(x = "Small cod density",
       y = "Flounder visregFit")

p2 <- ggplot(mfle_d2$fit, aes(x = large_cod_sc, y = visregFit)) +
  geom_point(aes(y = visregRes), data = mfle_d2$res, size = 1, alpha = 0.2, color = "grey30") +
  geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2) +
  geom_line(color = pal, size = 1.2) +
  labs(x = "Large cod density",
       y = "Flounder visregFit")

p3 <- ggplot(mscod_d$fit, aes(x = flounder_sc, y = visregFit)) +
  geom_point(aes(y = visregRes), data = mscod_d$res, size = 1, alpha = 0.2, color = "grey30") +
  geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2) +
  geom_line(color = pal, size = 1.2) +
  labs(x = "Flounder density",
       y = "Small cod visregFit")


p4 <- ggplot(mlcod_d$fit, aes(x = flounder_sc, y = visregFit)) +
  geom_point(aes(y = visregRes), data = mlcod_d$res, size = 1, alpha = 0.2, color = "grey30") +
  geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2) +
  geom_line(color = pal, size = 1.2) +
  labs(x = "Flounder density",
       y = "Large cod visregFit")

p1 + p2 + p3 + p4 & theme(aspect.ratio = 1)


ggsave("figures/density_coefs_smooth.pdf", width = 20, height = 20, units = "cm")

Plot coefficients

# Q1
# small cod
small_cod_q1 <- tidy(mcod_s_q1_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Small cod",
         quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# large cod
large_cod_q1 <- tidy(mcod_l_q1_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Large cod",
         quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# flounder
flounder_q1 <- tidy(mfle_q1_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Flounder",
         quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# Q4
# small cod
small_cod_q4 <- tidy(mcod_s_q4_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Small cod",
         quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# large cod
large_cod_q4 <- tidy(mcod_l_q4_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Large cod",
         quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# flounder
flounder_q4 <- tidy(mfle_q4_space, effects = "fixed", conf.int = TRUE) %>%
  mutate(species = "Flounder",
         quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA


coefs <- bind_rows(small_cod_q1,
                   large_cod_q1,
                   flounder_q1,
                      
                   small_cod_q4,
                   large_cod_q4,
                   flounder_q4) %>% 
  mutate(param_group = ifelse(term %in% c("fsubstratebedrock",
                                          "fsubstratehard clay",
                                          "fsubstratehard-bottom complex",
                                          "fsubstratemud",
                                          "fsubstratesand"),
                              "Substrate",
                              "Continious")) %>% 
  mutate(term = ifelse(term == "fsubstratebedrock", "Bedrock", term),
         term = ifelse(term == "fsubstratehard clay", "Hard clay", term),
         term = ifelse(term == "fsubstratehard-bottom complex", "Hard-bottom\ncomplex", term),
         term = ifelse(term == "fsubstratemud", "Mud", term),
         term = ifelse(term == "fsubstratesand", "Sand", term),
         term = ifelse(term == "depth_sc", "Depth", term),
         term = ifelse(term == "depth_sq_sc", "Depth^2", term),
         term = ifelse(term == "temp_sc", "Temperature", term),
         term = ifelse(term == "temp_sq_sc", "Temperature^2", term),
         term = ifelse(term == "oxy_sc", "Oxygen", term),
         term = ifelse(term == "sal_sc", "Salinity", term))
#> mutate: new variable 'param_group' (character) with 2 unique values and 0% NA
#> mutate: changed 66 values (29%) of 'term' (0 new NA)

coefs %>% 
  filter(!grepl('year', term)) %>% 
  ggplot(aes(term, estimate, color = species, shape = factor(quarter))) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  labs(shape = "Quarter") + 
  facet_wrap(~param_group, scales = "free") +
  labs(x = "Predictor", y = "Standardized coefficient") +
  theme_plot() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90),
        aspect.ratio = 1) +
  NULL 
#> filter: removed 162 rows (71%), 66 rows remaining


ggsave("figures/habitat_coefs.pdf", width = 20, height = 20, units = "cm")

Predict on grid

mcod_s_q1_pred <- predict(mcod_s_q1_space, newdata = pred_grid_q1)
mcod_s_q4_pred <- predict(mcod_s_q4_space, newdata = pred_grid_q4)

mcod_l_q1_pred <- predict(mcod_l_q1_space, newdata = pred_grid_q1)
mcod_l_q4_pred <- predict(mcod_l_q4_space, newdata = pred_grid_q4)

mfle_q1_pred <- predict(mfle_q1_space, newdata = pred_grid_q1)
mfle_q4_pred <- predict(mfle_q4_space, newdata = pred_grid_q4)

Plot predictions on map

p1 <- plot_map + 
  geom_raster(data = bind_rows(mcod_s_q1_pred, mcod_s_q4_pred) %>%
                filter(year %in% c("1995", "2019")),
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_grid(quarter~year) +
  ggtitle("Cod < 25 cm") +
  scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining

p2 <- plot_map + 
  geom_raster(data = bind_rows(mcod_l_q1_pred, mcod_l_q4_pred) %>%
                filter(year %in% c("1995", "2019")),
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_grid(quarter~year) +
  ggtitle("Cod > 25 cm") +
  scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining

p3 <- plot_map + 
  geom_raster(data = bind_rows(mfle_q1_pred, mfle_q4_pred) %>%
                filter(year %in% c("1995", "2019")),
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_grid(quarter~year) +
  ggtitle("Flounder") +
  scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining

((p1 + p2) / (p3 + plot_spacer())) + plot_annotation(tag_levels = "A") &
  theme(legend.text = element_text(angle = 90))


# Check flounder...

plot_map + 
  geom_raster(data = mfle_q1_pred,
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_wrap(~year) +
  scale_fill_viridis(trans = "sqrt")


plot_map + 
  geom_raster(data = mfle_q4_pred,
              aes(X*1000, Y*1000, fill = exp(est))) +
  facet_wrap(~year) +
  scale_fill_viridis(trans = "sqrt")

Calculate spatial overlap indices

all_pred_q1 <- mcod_s_q1_pred %>% mutate(small_cod = exp(est))
#> mutate: new variable 'small_cod' (double) with 296,380 unique values and 0% NA

all_pred_q1 <- all_pred_q1 %>%
  mutate(large_cod = exp(mcod_l_q1_pred$est),
         fle = exp(mfle_q1_pred$est))
#> mutate: new variable 'large_cod' (double) with 296,380 unique values and 0% NA
#>         new variable 'fle' (double) with 296,380 unique values and 0% NA

all_pred_q4 <- mcod_s_q4_pred %>% mutate(small_cod = exp(est))
#> mutate: new variable 'small_cod' (double) with 296,380 unique values and 0% NA

all_pred_q4 <- all_pred_q4 %>%
  mutate(large_cod = exp(mcod_l_q4_pred$est),
         fle = exp(mfle_q4_pred$est))
#> mutate: new variable 'large_cod' (double) with 296,380 unique values and 0% NA
#>         new variable 'fle' (double) with 296,380 unique values and 0% NA

all_pred <- bind_rows(all_pred_q1, all_pred_q4)

# Calculate biomass overlap
loc_collocspfn <- function(prey, pred) {
  p_prey <- prey/sum(prey, na.rm = T)
  p_pred <- pred/sum(pred, na.rm = T)
  (p_prey*p_pred)/(sqrt(sum(p_prey^2, na.rm = T))* sqrt(sum(p_pred^2, na.rm = T)))
}

loc_collocspfn_tot <- function(prey, pred) {
  p_prey <- prey/sum(prey, na.rm = T)
  p_pred <- pred/sum(pred, na.rm = T)
  sum((p_prey*p_pred))/(sqrt(sum(p_prey^2, na.rm = T))* sqrt(sum(p_pred^2, na.rm = T)))
}

all_pred <- all_pred %>% 
  group_by(year, quarter) %>% 
  mutate(scod_fle_ovr = loc_collocspfn(pred = small_cod, prey = fle),
         scod_fle_ovr_tot = loc_collocspfn_tot(pred = small_cod, prey = fle),
         lcod_fle_ovr = loc_collocspfn(pred = large_cod, prey = fle),
         lcod_fle_ovr_tot = loc_collocspfn_tot(pred = large_cod, prey = fle)) %>% 
  ungroup()
#> group_by: 2 grouping variables (year, quarter)
#> mutate (grouped): new variable 'scod_fle_ovr' (double) with 592,760 unique values and 0% NA
#>                   new variable 'scod_fle_ovr_tot' (double) with 56 unique values and 0% NA
#>                   new variable 'lcod_fle_ovr' (double) with 592,760 unique values and 0% NA
#>                   new variable 'lcod_fle_ovr_tot' (double) with 56 unique values and 0% NA
#> ungroup: no grouping variables

Plot

plot_map_fc +
  geom_raster(data = all_pred %>% filter(quarter == 1),
              aes(x = X, y = Y, fill = scod_fle_ovr)) +
  scale_fill_viridis_c(trans = "sqrt", name = "small cod-flounder") +
  facet_wrap(~ year, ncol = 5) +
  theme(legend.text = element_text(angle = 90)) + 
  NULL
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
#> Warning: Removed 296380 rows containing missing values (geom_raster).


plot_map_fc +
  geom_raster(data = all_pred %>% filter(quarter == 4),
              aes(x = X, y = Y, fill = lcod_fle_ovr)) +
  scale_fill_viridis_c(trans = "sqrt", name = "large cod-flounder") +
  facet_wrap(~ year, ncol = 5) +
  theme(legend.text = element_text(angle = 90)) + 
  NULL
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
#> Warning: Removed 296380 rows containing missing values (geom_raster).


# In space for selected year
all_pred2 <- all_pred %>%
  filter(year %in% c(1995, 2015)) %>% 
  dplyr::select(X, Y, scod_fle_ovr, lcod_fle_ovr, quarter, year) %>% 
  pivot_longer(c(scod_fle_ovr, lcod_fle_ovr)) %>% 
  rename(overlap = name) %>% 
  mutate(overlap = ifelse(overlap == "scod_fle_ovr", "Small cod : Flounder", overlap),
         overlap = ifelse(overlap == "lcod_fle_ovr", "Large cod : Flounder", overlap))
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
#> pivot_longer: reorganized (scod_fle_ovr, lcod_fle_ovr) into (name, value) [was 42340x6, now 84680x6]
#> rename: renamed one variable (overlap)
#> mutate: changed 84,680 values (100%) of 'overlap' (0 new NA)

# Since it's mainly an difference in elevation between quarters, I just do it for Q1 here
p1 <- plot_map_fc +
  geom_raster(data = filter(all_pred2, quarter == 1),
              aes(x = X*1000, y = Y*1000, fill = value)) +
  scale_fill_viridis_c(trans = "sqrt", name = "Overlap") +
  facet_grid(year ~ overlap) +
  theme(legend.text = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.25, 'cm'),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(size = 10)) + 
  ggtitle("Quarter 1") +
  NULL
#> filter: removed 42,340 rows (50%), 42,340 rows remaining

p2 <- plot_map_fc +
  geom_raster(data = filter(all_pred2, quarter == 4),
              aes(x = X*1000, y = Y*1000, fill = value)) +
  scale_fill_viridis_c(trans = "sqrt", name = "Overlap") +
  facet_grid(year ~ overlap) +
  theme(legend.text = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.25, 'cm'),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(size = 10)) +
  ggtitle("Quarter 4") +
  NULL
#> filter: removed 42,340 rows (50%), 42,340 rows remaining

# Over space
p5 <- all_pred %>% 
  distinct(year, quarter, .keep_all = TRUE) %>% 
  ggplot(aes(year, lcod_fle_ovr_tot, color = factor(quarter), fill = factor(quarter))) + 
  labs(color = "Quarter", x = "Year", y = "Overlap index") + 
  guides(fill = "none") +
  stat_smooth(method = "gam", formula = y~s(x, k = 3)) +
  geom_point() + 
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#> distinct: removed 592,704 rows (>99%), 56 rows remaining

p6 <- all_pred %>% 
  distinct(year, quarter, .keep_all = TRUE) %>% 
  ggplot(aes(year, scod_fle_ovr_tot, color = factor(quarter), fill = factor(quarter))) + 
  labs(color = "Quarter", x = "Year", y = "Overlap index") + 
  guides(fill = "none") +
  stat_smooth(method = "gam", formula = y~s(x, k = 3)) +
  geom_point() +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#> distinct: removed 592,704 rows (>99%), 56 rows remaining

((p1 | p2) / ((p5 | p6) + plot_layout(guides = "collect"))) + plot_layout(heights = c(1.1, 1)) + plot_annotation(tag_levels = "A") & theme(aspect.ratio = 1)


ggsave("figures/spatial_overlap.pdf", width = 20, height = 20, units = "cm")

Conditional effects

# Scale with respect to data!
nd_depth <- data.frame(depth = seq(min(cpue2$depth), max(cpue2$depth), length.out = 100)) %>% 
  mutate(depth_ct = depth - mean(cpue2$depth),
         depth_sq = depth_ct * depth_ct,
         depth_sq_sc = (depth_sq - mean(cpue2$depth_sq)) / sd(cpue2$depth_sq),
         depth_sc = (depth - mean(cpue2$depth)) / sd(cpue2$depth),
         temp_sq_sc = 0,
         temp_sc = 0,
         oxy_sc = 0,
         sal_sc = 0,
         year = 1999L,
         fyear = as.factor(1999),
         fsubstrate = "mud")
#> mutate: new variable 'depth_ct' (double) with 100 unique values and 0% NA
#>         new variable 'depth_sq' (double) with 100 unique values and 0% NA
#>         new variable 'depth_sq_sc' (double) with 100 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 100 unique values and 0% NA
#>         new variable 'temp_sq_sc' (double) with one unique value and 0% NA
#>         new variable 'temp_sc' (double) with one unique value and 0% NA
#>         new variable 'oxy_sc' (double) with one unique value and 0% NA
#>         new variable 'sal_sc' (double) with one unique value and 0% NA
#>         new variable 'year' (integer) with one unique value and 0% NA
#>         new variable 'fyear' (factor) with one unique value and 0% NA
#>         new variable 'fsubstrate' (character) with one unique value and 0% NA

# Q1
margin_small_cod_pred_q1_depth <- predict(mcod_s_q1_space,
                                          newdata = nd_depth,
                                          se_fit = TRUE,
                                          re_form = NA) %>%
  mutate(species = "small_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_large_cod_pred_q1_depth <- predict(mcod_l_q1_space,
                                          newdata = nd_depth,
                                          se_fit = TRUE,
                                          re_form = NA) %>%
  mutate(species ="large_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_fle_pred_q1_depth <- predict(mfle_q1_space,
                                    newdata = nd_depth,
                                    se_fit = TRUE,
                                    re_form = NA) %>%
  mutate(species = "flounder", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# Q4
margin_small_cod_pred_q4_depth <- predict(mcod_s_q4_space,
                                          newdata = nd_depth,
                                          se_fit = TRUE,
                                          re_form = NA) %>%
  mutate(species = "small_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_large_cod_pred_q4_depth <- predict(mcod_l_q4_space,
                                          newdata = nd_depth,
                                          se_fit = TRUE,
                                          re_form = NA) %>%
  mutate(species = "large_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_fle_pred_q4_depth <- predict(mfle_q4_space,
                                    newdata = nd_depth,
                                    se_fit = TRUE,
                                    re_form = NA) %>%
  mutate(species = "flounder", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margins_depth <- bind_rows(margin_small_cod_pred_q1_depth,
                           margin_large_cod_pred_q1_depth,
                           margin_fle_pred_q1_depth,
                           margin_small_cod_pred_q4_depth,
                           margin_large_cod_pred_q4_depth,
                           margin_fle_pred_q4_depth) %>% 
    mutate(species = ifelse(species == "flounder", "Flounder", species),
           species = ifelse(species == "small_cod", "Small cod", species),
           species = ifelse(species == "large_cod", "Large cod", species))
#> mutate: changed 600 values (100%) of 'species' (0 new NA)

# Temperature
nd_temp <- data.frame(temp = seq(min(cpue2$temp), max(cpue2$temp), length.out = 100)) %>% 
  mutate(temp_ct = temp - mean(cpue2$temp),
         temp_sq = temp_ct * temp_ct,
         temp_sq_sc = (temp_sq - mean(cpue2$temp_sq)) / sd(cpue2$temp_sq),
         temp_sc = (temp - mean(cpue2$temp)) / sd(cpue2$temp),
         depth_sq_sc = 0,
         depth_sc = 0,
         oxy_sc = 0,
         sal_sc = 0,
         year = 1999L,
         fyear = as.factor(1999),
         fsubstrate = "mud")
#> mutate: new variable 'temp_ct' (double) with 100 unique values and 0% NA
#>         new variable 'temp_sq' (double) with 100 unique values and 0% NA
#>         new variable 'temp_sq_sc' (double) with 100 unique values and 0% NA
#>         new variable 'temp_sc' (double) with 100 unique values and 0% NA
#>         new variable 'depth_sq_sc' (double) with one unique value and 0% NA
#>         new variable 'depth_sc' (double) with one unique value and 0% NA
#>         new variable 'oxy_sc' (double) with one unique value and 0% NA
#>         new variable 'sal_sc' (double) with one unique value and 0% NA
#>         new variable 'year' (integer) with one unique value and 0% NA
#>         new variable 'fyear' (factor) with one unique value and 0% NA
#>         new variable 'fsubstrate' (character) with one unique value and 0% NA

# Q1
margin_small_cod_pred_q1_temp <- predict(mcod_s_q1_space,
                                         newdata = nd_temp,
                                         se_fit = TRUE,
                                         re_form = NA) %>%
  mutate(species = "small_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_large_cod_pred_q1_temp <- predict(mcod_l_q1_space,
                                         newdata = nd_temp,
                                         se_fit = TRUE,
                                         re_form = NA) %>%
  mutate(species ="large_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_fle_pred_q1_temp <- predict(mfle_q1_space,
                                   newdata = nd_temp,
                                   se_fit = TRUE,
                                   re_form = NA) %>%
  mutate(species = "flounder", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

# Q4
margin_small_cod_pred_q4_temp <- predict(mcod_s_q4_space,
                                         newdata = nd_temp,
                                         se_fit = TRUE,
                                         re_form = NA) %>%
  mutate(species = "small_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_large_cod_pred_q4_temp <- predict(mcod_l_q4_space,
                                         newdata = nd_temp,
                                         se_fit = TRUE,
                                         re_form = NA) %>%
  mutate(species = "large_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margin_fle_pred_q4_temp <- predict(mfle_q4_space,
                                   newdata = nd_temp,
                                   se_fit = TRUE,
                                   re_form = NA) %>%
  mutate(species = "flounder", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#>         new variable 'quarter' (double) with one unique value and 0% NA

margins_temp <- bind_rows(margin_small_cod_pred_q1_temp,
                          margin_large_cod_pred_q1_temp,
                          margin_fle_pred_q1_temp,
                          margin_small_cod_pred_q4_temp,
                          margin_large_cod_pred_q4_temp,
                          margin_fle_pred_q4_temp) %>% 
    mutate(species = ifelse(species == "flounder", "Flounder", species),
           species = ifelse(species == "small_cod", "Small cod", species),
           species = ifelse(species == "large_cod", "Large cod", species))
#> mutate: changed 600 values (100%) of 'species' (0 new NA)

# Plot!
p1 <- ggplot(margins_depth, 
             aes(depth, exp(est), ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se),
                 color = species, fill = species)) +
  geom_line() +
  facet_wrap(~quarter, scales = "free") +
  geom_ribbon(alpha = 0.4, color = NA) +
  coord_cartesian(xlim = c(10, 170)) +
  theme(aspect.ratio = 1) +
  labs(x = "Depth (m)", y = "Biomass density", color = "", fill = "") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)


p2 <- ggplot(margins_temp,
             aes(temp, exp(est), ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se),
                 color = species, fill = species)) +
  geom_line() +
  facet_wrap(~quarter, scales = "free") +
  geom_ribbon(alpha = 0.4, color = NA) +
  coord_cartesian(xlim = c(1, 14)) +
  theme(aspect.ratio = 1) +
  labs(x = "Temperature (°C)", y = "Biomass density", color = "", fill = "") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)

(p1 / p2) + plot_annotation(tag_levels = "A") + plot_layout(guides = "collect")


# ggsave("figures/supp/conditional_effects.pdf", width = 20, height = 20, units = "cm#")

# Scaled
p3 <- margins_temp %>% 
  ungroup() %>% 
  group_by(species, quarter) %>% 
  mutate(max_est = max(exp(est)),
         est_sc = exp(est) / max_est) %>% 
  ggplot(aes(temp, est_sc, color = species, fill = species)) +
  geom_line(size = 1.3) +
  facet_wrap(~quarter, scales = "free") +
  coord_cartesian(xlim = c(1, 14)) +
  theme(aspect.ratio = 1) +
  labs(x = "Temperature (°C)", y = "Scaled biomass density", color = "", fill = "") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> ungroup: no grouping variables
#> group_by: 2 grouping variables (species, quarter)
#> mutate (grouped): new variable 'max_est' (double) with 6 unique values and 0% NA
#>                   new variable 'est_sc' (double) with 595 unique values and 0% NA

p4 <- margins_depth %>% 
  ungroup() %>% 
  group_by(species, quarter) %>% 
  mutate(max_est = max(exp(est)),
         est_sc = exp(est) / max_est) %>% 
  ggplot(aes(depth, est_sc, color = species, fill = species)) +
  geom_line(size = 1.3) +
  facet_wrap(~quarter, scales = "free") +
  coord_cartesian(xlim = c(10, 170)) +
  theme(aspect.ratio = 1) +
  labs(x = "Depth (m)", y = "Scaled biomass density", color = "", fill = "") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> ungroup: no grouping variables
#> group_by: 2 grouping variables (species, quarter)
#> mutate (grouped): new variable 'max_est' (double) with 6 unique values and 0% NA
#>                   new variable 'est_sc' (double) with 595 unique values and 0% NA

(p3 / p4) + plot_annotation(tag_levels = "A") + plot_layout(guides = "collect")


# ggsave("figures/conditional_effects.pdf", width = 20, height = 20, units = "cm")

Weighted quantiles

head(all_pred) %>% as.data.frame()
#>     X    Y depth year      lon      lat substrate quarter      oxy     temp
#> 1 450 5984    11 1993 14.23718 54.00188      sand       1 8.713132 2.821855
#> 2 454 5984    11 1993 14.29820 54.00225      sand       1 8.743857 2.811376
#> 3 458 5984    11 1993 14.35922 54.00259      sand       1 8.745697 2.806495
#> 4 462 5984    12 1993 14.42024 54.00290      sand       1 8.759304 2.802808
#> 5 466 5984    12 1993 14.48127 54.00318      sand       1 8.765212 2.805644
#> 6 470 5984    11 1993 14.54229 54.00343      sand       1 8.805174 2.792613
#>        sal sub_div ices_rect density_saduria biomass_spr biomass_her
#> 1 7.293935      24      37G4               0          NA          NA
#> 2 7.346682      24      37G4               0          NA          NA
#> 3 7.383028      24      37G4               0          NA          NA
#> 4 7.418652      24      37G4               0          NA          NA
#> 5 7.501428      24      37G4               0          NA          NA
#> 6 7.585449      24      37G4               0          NA          NA
#>   biomass_spr_sd biomass_her_sd  depth_ct depth_sq depth_sq_sc  depth_sc
#> 1             NA             NA -36.63035 1341.783   0.4698617 -1.343238
#> 2             NA             NA -36.63035 1341.783   0.4698617 -1.343238
#> 3             NA             NA -36.63035 1341.783   0.4698617 -1.343238
#> 4             NA             NA -35.63035 1269.522   0.4131041 -1.306567
#> 5             NA             NA -35.63035 1269.522   0.4131041 -1.306567
#> 6             NA             NA -36.63035 1341.783   0.4698617 -1.343238
#>     temp_ct  temp_sq temp_sq_sc   temp_sc   oxy_sc    sal_sc fyear fsubstrate
#> 1 -3.484890 12.14446  0.5884181 -1.305532 1.178372 -1.098960  1993       sand
#> 2 -3.495369 12.21760  0.5969916 -1.309458 1.190285 -1.082916  1993       sand
#> 3 -3.500250 12.25175  0.6009941 -1.311287 1.190999 -1.071861  1993       sand
#> 4 -3.503937 12.27757  0.6040213 -1.312668 1.196275 -1.061026  1993       sand
#> 5 -3.501101 12.25771  0.6016927 -1.311606 1.198566 -1.035849  1993       sand
#> 6 -3.514132 12.34912  0.6124078 -1.316487 1.214061 -1.010293  1993       sand
#>         est est_non_rf     est_rf    omega_s epsilon_st  small_cod large_cod
#> 1 -3.128299  -2.913320 -0.2149794 -0.3565650  0.1415856 0.04379222  6.078414
#> 2 -3.475978  -2.915583 -0.5603953 -0.7218136  0.1614182 0.03093157  5.816645
#> 3 -3.803139  -2.919747 -0.8833924 -1.0647534  0.1813610 0.02230065  5.483502
#> 4 -3.886795  -2.811876 -1.0749196 -1.2701738  0.1952542 0.02051097  5.240582
#> 5 -4.076864  -2.810417 -1.2664467 -1.4755941  0.2091473 0.01696057  4.760898
#> 6 -4.302368  -2.921767 -1.3806012 -1.6096988  0.2290976 0.01353646  4.322854
#>        fle scod_fle_ovr scod_fle_ovr_tot lcod_fle_ovr lcod_fle_ovr_tot
#> 1 13.58323 6.025629e-08        0.3323209 9.576900e-07        0.5491266
#> 2 13.07450 4.096655e-08        0.3323209 8.821231e-07        0.5491266
#> 3 12.76604 2.883873e-08        0.3323209 8.119809e-07        0.5491266
#> 4 13.48553 2.801926e-08        0.3323209 8.197456e-07        0.5491266
#> 5 14.17812 2.435911e-08        0.3323209 7.829590e-07        0.5491266
#> 6 14.43847 1.979833e-08        0.3323209 7.239743e-07        0.5491266

all_pred_sub <- all_pred %>%
  dplyr::select(small_cod, large_cod, fle, year, temp, oxy, sal, depth, quarter) %>% 
  pivot_longer(c(small_cod, large_cod, fle)) %>% 
  rename(Species = name, density = value) 
#> pivot_longer: reorganized (small_cod, large_cod, fle) into (name, value) [was 592760x9, now 1778280x8]
#> rename: renamed 2 variables (Species, density)

wm <- all_pred_sub %>%
  pivot_longer(c(temp, oxy, sal, depth)) %>% 
  mutate(name = ifelse(name == "temp", "Temperature (°C)", name),
         name = ifelse(name == "sal", "Salinity", name),
         name = ifelse(name == "depth", "Depth", name),
         name = ifelse(name == "oxy", "Oxygen (ml/L)", name)) %>% 
  group_by(year, quarter, Species, name) %>%
  summarise("Weighted median" = hutils::weighted_quantile(v = value, w = density, p = c(0.5)),
            "1st quartile" = hutils::weighted_quantile(v = value, w = density, p = c(0.25)),
            "3rd quartile" = hutils::weighted_quantile(v = value, w = density, p = c(0.75))) %>% 
  rename(env_var = name) %>% 
  mutate(Species = ifelse(Species == "fle", "Flounder", Species),
         Species = ifelse(Species == "large_cod", "Large cod", Species),
         Species = ifelse(Species == "small_cod", "Small cod", Species))
#> pivot_longer: reorganized (temp, oxy, sal, depth) into (name, value) [was 1778280x8, now 7113120x6]
#> mutate: changed 7,113,120 values (100%) of 'name' (0 new NA)
#> group_by: 4 grouping variables (year, quarter, Species, name)
#> summarise: now 672 rows and 7 columns, 3 group variables remaining (year, quarter, Species)
#> rename: renamed one variable (env_var)
#> mutate (grouped): changed 672 values (100%) of 'Species' (0 new NA)

ggplot(wm, aes(year, `Weighted median`, color = Species, fill = Species)) +
  geom_line() + 
  ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)


ggplot(wm, aes(year, `Weighted median`, color = Species)) +
  geom_line(size = 0.8) + 
  labs(linetype = "Quarter", x = "Year") +
  ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)


#ggsave("figures/weighted_quantiles.pdf", width = 20, height = 20, units = "cm")

# Scaled
wm %>% 
  group_by(Species, quarter, env_var) %>% 
  mutate(`Weighted median scaled` = `Weighted median` - mean(`Weighted median`)) %>% 
  ggplot(aes(year, `Weighted median scaled`, color = Species)) +
  geom_line(size = 0.8) + 
  labs(linetype = "Quarter", x = "Year") +
  ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
  scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
  scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> group_by: 3 grouping variables (Species, quarter, env_var)
#> mutate (grouped): new variable 'Weighted median scaled' (double) with 656 unique values and 0% NA


#ggsave("figures/supp/weighted_quantiles_scaled.pdf", width = 20, height = 20, units = "cm")

# ggplot(wm, aes(year, Median, color = Species, fill = Species)) +
#   geom_ribbon(aes(ymin = `1st quartile`, ymax = `3rd quartile`), alpha = 0.4, color = NA) + 
#   facet_wrap(env_var ~ quarter, scales = "free", ncol = 2) + 
#   scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
#   scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
# 
# ggsave("figures/supp/weighted_quantiles_range.pdf", width = 20, height = 20, units = "cm")

[Extra for bentfish - calculate biomass index per subdivision]

Compare with Q4 from previous index

cpue3 <- cpue2 %>%
  mutate(tot_cod = large_cod + small_cod) %>%
  filter(quarter == 4)
#> mutate: new variable 'tot_cod' (double) with 7,870 unique values and 0% NA
#> filter: removed 5,271 rows (59%), 3,595 rows remaining

spde <- make_mesh(cpue3, xy_cols = c("X", "Y"),
                  n_knots = 200, 
                  type = "kmeans", seed = 42)

cod_ind <- sdmTMB(tot_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                    temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                  data = cpue3,
                  mesh = spde,
                  family = tweedie(link = "log"),
                  spatiotemporal = "AR1",
                  spatial = "on",
                  time = "year",
                  reml = FALSE,
                  control = sdmTMBcontrol(newton_loops = 1))

fle_ind <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
                    temp_sc + temp_sq_sc + oxy_sc + sal_sc,
                  data = cpue3,
                  mesh = spde,
                  family = tweedie(link = "log"),
                  spatiotemporal = "AR1",
                  spatial = "on",
                  time = "year",
                  reml = FALSE,
                  control = sdmTMBcontrol(newton_loops = 1))

sanity(cod_ind)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
sanity(fle_ind)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
# cod_ind2 <- sdmTMB(tot_cod ~ 0 + fyear + s(depth_sc),
#                   data = cpue3,
#                   mesh = spde,
#                   family = tweedie(link = "log"),
#                   spatiotemporal = "AR1",
#                   spatial = "on",
#                   time = "year",
#                   reml = FALSE,
#                   control = sdmTMBcontrol(newton_loops = 1))
# 
# fle_ind2 <- sdmTMB(flounder ~ 0 + fyear + s(depth_sc),
#                   data = cpue3,
#                   mesh = spde,
#                   family = tweedie(link = "log"),
#                   spatiotemporal = "AR1",
#                   spatial = "on",
#                   time = "year",
#                   reml = FALSE,
#                   control = sdmTMBcontrol(newton_loops = 1))
# 
# sanity(cod_ind2)
# sanity(fle_ind2)

Now try old data as well

# d_old <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue_q_1_4.csv")
# 
# # Filter to match the data I want to predict on
# d_old <- d_old %>%
#   filter(year > 1992 & year < 2020 & quarter == 4)
# 
# # Calculate standardized variables
# d_old <- d_old %>% 
#   mutate(depth_sc = (depth - mean(depth))/sd(depth),
#          X = X/1000,
#          Y = Y/1000,
#          year = as.integer(year),
#          fyear  = as.factor(year)) %>% 
#   drop_na(depth) %>% 
#   rename("density_cod" = "density") # to fit better with how flounder is named
# 
# nrow(d_old)
# nrow(cpue3)
# 
# # Plot comparison
# p1 <- ggplot() + 
#   geom_histogram(data = d_old, aes(density_cod, fill = "old"), alpha = 0.5) +
#   geom_histogram(data = cpue3, aes(tot_cod, fill = "new"), alpha = 0.5) + 
#   scale_y_continuous(trans = "log")
# 
# p2 <- ggplot() + 
#   geom_histogram(data = d_old, aes(density_fle, fill = "old"), alpha = 0.5) +
#   geom_histogram(data = cpue3, aes(flounder, fill = "new"), alpha = 0.5) + 
#   scale_y_continuous(trans = "log")
# 
# p1 + p2
# spde_old <- make_mesh(d_old, xy_cols = c("X", "Y"),
#                   n_knots = 200, 
#                   type = "kmeans", seed = 42)
# 
# cod_ind3 <- sdmTMB(density_cod ~ 0 + fyear + s(depth_sc),
#                   data = d_old,
#                   mesh = spde_old,
#                   family = tweedie(link = "log"),
#                   spatiotemporal = "AR1",
#                   spatial = "on",
#                   time = "year",
#                   reml = FALSE,
#                   control = sdmTMBcontrol(newton_loops = 1))
# 
# fle_ind3 <- sdmTMB(density_fle ~ 0 + fyear + s(depth_sc),
#                   data = d_old,
#                   mesh = spde_old,
#                   family = tweedie(link = "log"),
#                   spatiotemporal = "AR1",
#                   spatial = "on",
#                   time = "year",
#                   reml = FALSE,
#                   control = sdmTMBcontrol(newton_loops = 1))
# 
# sanity(cod_ind3)
# sanity(fle_ind3)
### Cod
# SD 24
preds_cod24_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 513,744 rows (87%), 79,016 rows remaining

# SD 25
preds_cod25_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 446,040 rows (75%), 146,720 rows remaining

# SD 26
preds_cod26_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 457,856 rows (77%), 134,904 rows remaining

# SD 27
preds_cod27_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 498,120 rows (84%), 94,640 rows remaining

# SD 28
preds_cod28_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 455,280 rows (77%), 137,480 rows remaining

# Now calculate the CPUE index (average)
index24_cod_sim <- get_index_sims(preds_cod24_sim, area = rep(4*4, nrow(preds_cod24_sim)))
index25_cod_sim <- get_index_sims(preds_cod25_sim, area = rep(4*4, nrow(preds_cod25_sim)))
index26_cod_sim <- get_index_sims(preds_cod26_sim, area = rep(4*4, nrow(preds_cod26_sim)))
index27_cod_sim <- get_index_sims(preds_cod27_sim, area = rep(4*4, nrow(preds_cod27_sim)))
index28_cod_sim <- get_index_sims(preds_cod28_sim, area = rep(4*4, nrow(preds_cod28_sim)))

# Add some columns
index24_cod_sim <- index24_cod_sim %>% mutate(sub_div = "24", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index25_cod_sim <- index25_cod_sim %>% mutate(sub_div = "25", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index26_cod_sim <- index26_cod_sim %>% mutate(sub_div = "26", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index27_cod_sim <- index27_cod_sim %>% mutate(sub_div = "27", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index28_cod_sim <- index28_cod_sim %>% mutate(sub_div = "28", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA


### Flounder
# SD 24
preds_fle24_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 513,744 rows (87%), 79,016 rows remaining

# SD 25
preds_fle25_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 446,040 rows (75%), 146,720 rows remaining

# SD 26
preds_fle26_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 457,856 rows (77%), 134,904 rows remaining

# SD 27
preds_fle27_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 498,120 rows (84%), 94,640 rows remaining

# SD 28
preds_fle28_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 455,280 rows (77%), 137,480 rows remaining

# Now calculate the CPUE index (average)
index24_fle_sim <- get_index_sims(preds_fle24_sim, area = rep(4*4, nrow(preds_fle24_sim)))
index25_fle_sim <- get_index_sims(preds_fle25_sim, area = rep(4*4, nrow(preds_fle25_sim)))
index26_fle_sim <- get_index_sims(preds_fle26_sim, area = rep(4*4, nrow(preds_fle26_sim)))
index27_fle_sim <- get_index_sims(preds_fle27_sim, area = rep(4*4, nrow(preds_fle27_sim)))
index28_fle_sim <- get_index_sims(preds_fle28_sim, area = rep(4*4, nrow(preds_fle28_sim)))

# Add some columns
index24_fle_sim <- index24_fle_sim %>% mutate(sub_div = "24", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index25_fle_sim <- index25_fle_sim %>% mutate(sub_div = "25", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index26_fle_sim <- index26_fle_sim %>% mutate(sub_div = "26", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index27_fle_sim <- index27_fle_sim %>% mutate(sub_div = "27", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index28_fle_sim <- index28_fle_sim %>% mutate(sub_div = "28", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

div_index_sim <- bind_rows(index24_cod_sim,
                           index25_cod_sim,
                           index26_cod_sim,
                           index27_cod_sim,
                           index28_cod_sim,
                           index24_fle_sim,
                           index25_fle_sim,
                           index26_fle_sim,
                           index27_fle_sim,
                           index28_fle_sim) %>% 
  mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000) # convert to tonnes 
#> mutate: new variable 'est_t' (double) with 280 unique values and 0% NA
#>         new variable 'lwr_t' (double) with 280 unique values and 0% NA
#>         new variable 'upr_t' (double) with 280 unique values and 0% NA

Plot the index and save as csv

# # First compare with previous index:
# cod_fle_index_old <- read.csv("/Users/maxlindmark/Dropbox/Max work/R/cod_condition/output/cod_fle_index.csv") %>% 
#   mutate(sub_div = as.character(sub_div)) %>% 
#   filter(sub_div %in% c("25", "28"))
# 
# index_comp <- bind_rows(cod_fle_index_old %>% mutate(source = "cod_condition"),
#                         div_index_sim %>% mutate(source = "cod_interactions"))
# 
# ggplot(index_comp, aes(year, est_t, color = source, fill = source)) +
#   facet_wrap(sub_div ~ species, scales = "free") +
#   geom_line() +
#   geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t), alpha = 0.2, color = NA)
# 
# 
# ## If it isn't the same with the data, plot pred grids against each other

ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
  geom_line() +
  facet_wrap(~species, scales = "free") +
  geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2, color = NA) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL


div_index_sim %>% 
  filter(sub_div %in% c(25, 28)) %>% 
  ggplot(aes(year, est_t/1000, color = species, fill = species)) +
  geom_line() +
  facet_wrap(~sub_div, scales = "free", ncol = 1) +
  geom_ribbon(aes(year, ymin = lwr_t/1000, ymax = upr_t/1000), alpha = 0.2, color = NA) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "Species") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL
#> filter: removed 168 rows (60%), 112 rows remaining


ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
  geom_line() +
  facet_wrap(~species, scales = "free") +
  #geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL


ggplot(div_index_sim, aes(year, est_t, color = species)) +
  geom_line() +
  facet_wrap(~sub_div, scales = "free") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL


write.csv(div_index_sim, "output/cod_fle_index.csv")